Correlation-induced triplet superconductivity on the graphene lattice 
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We investigate the possibility of superconductivity on the graphene lattice within the repulsive 
Hubbard model using the variational cluster approximation (VGA). We find that singlet supercon- 
ductivity is impossible; instead, triplet superconductivity is favored, with four solutions that are 
close to each other in energy and differ by their symmetry; one is f-wave, the other three p-wave, 
including a p+ip solution that breaks time-reversal invariance. 



The electronic properties of graphene have been the 
object of much research since isolated sheets have been 
manipulated and characterized in 2004 (for a review, see 
[1]). Graphene has a very high mobility, displays an 
anomalous quantum Hall effect (at very high field) even 
at room temperature and also shows a universal conduc- 
tivity. The possibility of doping graphene by applying an 
electric field offers the prospect of carbon-based electron- 
ics. From a theoretical point of view, a peculiar feature 
of graphene is the cone-like dispersion relation around 
two unequivalent points in the Brillouin zone, which al- 
lows a low-energy description in terms of fermions obey- 
ing the (2 -|- 1) -dimensional Dirac equation. There have 
been speculations that a modified graphene system, ob- 
tained for instance by stacking planes or by disorder, 
could display magnetic or even superconducing order at 
room temperature [2-5]. The goal of this work is to check 
whether electronic correlations alone, as described by a 
repulsive Hubbard model, can induce superconductivity 
on doped graphene; thus we will ignore electron-phonon 
and long-range Coulomb interactions. Neglecting the lat- 
ter is likely a poor approximation very close to half-filling, 
since screening is then hindered by a small density of 
states, and we expect our work to be relevant mostly 
away from that point. 

Our study focuses on the two-dimensional Hubbard 




FIG. 1: (Golor online) 6- and fO-site clusters used in this 
work. The superlattice vectors are indicated by red arrows. 
The A and B sublattices are indicated by black and blue dots, 
respectively. 



model defined on the honeycomb lattice with nearest- 
neighbor hopping amplitude t, representing interac- 
tions within the tt band: In pure graphene and three- 
dimensional graphite, this band is half-filled, t is esti- 
mated at 2.8 eV and the on-site Coulomb repulsion U 
is expected to be around 17 eV [1]. We will set U = 6t 
in what follows. The honeycomb lattice is bipartite: it 
can be viewed as the sum of two interspersed triangular 
sublattices (A and B, see Fig. 1). The Hubbard model 
Hamiltonian is expressed as 
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where ci]a destroys (creates) and electron of spin a in 

a Wannier orbital at site r, and rir^a = ci]l-Cr^a is the 
number of electrons of spin a at site r. The three vectors 
ei,2.3 link a site of sublattice A with its three nearest 
neighbors (NN) on sublattice B, and are oriented at 120° 
of each other. The first sum runs over sites of the A sub- 
lattice only and contains all hopping terms. The second 
sum (the local Coulomb repulsion) runs over all sites. 

Superconductivity is characterized by the condensa- 
tion of Cooper pairs, which translates into a nonzero 
expectation value for some pairing operator. In mo- 
mentum space, it is convenient to arrange the cre- 
ation and annihilation operators for the two sublat- 
tices and the two spins into a four-component object 

/ = (cA,k,T>cs,k,T,c^ -k.i'C^ _k,i)- Pairing may occur 
either in the singlet or triplet channel; we will show be- 
low that the latter is favored. In the mean-field approx- 
imation, the Hamiltonian of the system would take the 
form Hmf = P'Hf . In the triplet channel, Ti. is 
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where /i is the chemical potential, — t J2j=i 2 3 e*'^ ''^ 
is the hopping function and 6'k and 77k describe pairing 
amphtudes between electrons of the same sublattice and 
different sublattices, respectively (this parametrization 
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assumes the pairings and hopping terms to be real when 
expressed in real space). 

It is a straightforward matter to show that the (four- 



Owing to the two bands of graphene, this expression does 
not have the classic BCS shape with an easily identifiable 
gap function Ak (the same goes for the singlet channel 
dispersion). The question of gap symmetry may instead 
be investigated by considering the dominant pairings in 
real space (see also Fig. 3 below) . The repulsive local in- 
teraction excludes the possibility of on-site pairing (i.e., 
has no constant term). Pairing operators may be de- 
fined between nearest-neighbor (NN) sites {i = 1,2,3), 
either in the singlet (— sign) or triplet (-|- sign) channel; 

The associated pairing functions are simply 77k = 
Sj=i 2 3*^3^^^^' = 0, where we introduced pair- 

ing amplitudes ai,2,3 in the three NN directions. In 
Ref. [2], the spin-singlet NN pairing Si was referred to 
as 'p-wave'; we will refrain from using this terminology 
since it might be confusing in the context of a point-group 
symmetry description of the possible pairing states. 

The honeycomb lattice is characterized by a Cqv sym- 
metry about the center of the hexag ons (^Cqv ^ Dq in 
two dimensions). A straightforward analysis in terms of 
group projection operators shows that the six pairings op- 
erators (4) can be arranged into four different irreducible 
representations of Cq^: 

Aa, = a{Si +S2 + S3) 

Ab, = Q:(ri +T2 + T3) 

Ae, = a{T, - T2) + P{T2 - T3) ^' 
Ae, = a{Si - S2) + P{S2 - S3) 

(the El and E2 representations are two-dimensional, and 
a, (3 are constants; the site index r suppressed). At the 
critical temperature, it is expected from Landau theory 
that the actual superconducting state fall into one of the 
above group representations. As the temperature is low- 
ered, additional symmetry breaking transitions may oc- 
cur so that the zero-temperature state may be a mixture 
of the above states (see also the symmetry analysis of 
Ref. [6] in the context of Cobaltates). 

The mean-field picture is useful to develop a physi- 
cal sense of the superconducting state, but we will not 



branch) dispersion relation derived from this mean-field 
Hamiltonian is 
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use it in computations. Instead, we use the more pow- 
erful variational cluster approximation (VGA), at zero- 
temperature. The VGA [7] is a variational method based 
on the electron self-energy, as defined in Potthoff's self- 
energy functional approach (SFA) [8] . The basic idea be- 
hind the SFA is to introduce a reference Hamiltonian 
H', with the same two-body interaction as the original 
Hamiltonian H, but with a different one-body part, so 
that H' may be solved exactly (numerically). In the 
VGA, the original lattice is tiled into a superlattice of 
identical clusters and H' is the cluster Hamiltonian; it 
differs from the original Hamiltonian H by the suppres- 
sion of intercluster hopping terms and the addition of 
Weiss fields associated with the broken symmetry phases 
of interest. The values of the Weiss fields, as well as other 
one-body parameters like the cluster's chemical potential, 
serve as variational parameters to optimize a functional 
fl whose expression is 

fi(t')=f]'-/ ^^lndet(l + (Go'-G^,-^)GO (6) 
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where Go is the non-interacting Green function of the 
original Hamiltonian, Gq the non-interacting Green func- 
tion of the cluster Hamiltonian and G' the exact Green 
function of the cluster Hamiltonian. At the physical self- 
energy, this functional approximates the grand potential 
Q. of the system. The only approximation comes from 
the limited space of self-energies on which the variational 
principle is applied, limited by the cluster size and by the 
number of variational parameters used. 

The VGA has been applied, for instance, to the prob- 
lem of competing phases in the high-Tc cuprates [9, 10] 
and in the layered organic conductors [11]. The VGA 
does not require any factorization of the interaction, and 
short-range correlations (within a cluster) are taken into 
account exactly. The Green function obtained from VGA 
is still defined on the infinite lattice. 

Fig. 1 illustrates the two clusters (L = 6 and L = 10 
sites, respectively) that were used in this work. While it 
is not possible in this case to perform any kind of finite- 
size scaling analysis, comparing the results from different 
clusters is useful in order to assess their robustness. The 
6-site cluster has the advantage of posessing the Gqu sym- 
metry of the lattice. The variational parameters used in 
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FIG. 2: (Color online) Sample Potthoff functional 57 as a 
function of Weiss field, for three triplet solutions and the cor- 
responding singlet solutions, at half-filling. The minima are 
shown as dots. Solution (9) is not shown since it involves two 
variational parameters. 



is this work are the (complex) coefhcients of the six dif- 
ferent pairing operators (4), as well as the cluster's chem- 
ical potential fj,' (including the latter in the variational 
set garantees thermodynamic consistency[10]). 

Four different superconducting solutions were found, 
all of them spin triplets. Thus the first conclusion of this 
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FIG. 3: (Color online) Density plots of the order parame- 
ter {cA,k,TCA,-k,i) as a function of k, for mean-field solutions 
having the same symmetry as the four triplet solutions found. 
The Brillouin zone (hexagon) is indicated. Negative (positive) 
values are represented by shades of blue (red); white is zero. 
Note how the topology of nodes (white lines) differs among 
the four solutions. The bottom-right panel only shows the 
real part of a complex order parameter; the imaginary part 
has conjugate nodes, so that the modulus of the order param- 
eter is constant along the Fermi surface. 



work is that singlet superconductivity does not occur in 
this system through a purely repulsive interaction. Each 
of the four solutions can be described by the coefficients 
(q!i, a2, CK3) of the triplet paring operator aiTi + 0:212 + 

{a, a, a) Si representation (f) (7) 

(a, —a, 0) El representation (p) (8) 

{a,a,—P) Si — i?i mixture (p) (9) 

(a,ia,0) broken T-reversal (p + ip) (10) 

Solutions (8), (9) and (10) also exist in rotated form, ob- 
tained for instance by permuting the coefficients. Fig. 2 
illustrates the dependence of the Potthoff functional on 
the Weiss field a for three of the above solutions, as well 
as for the corresponding trial singlet-SC states. The only 
extrema of the singlet states occurs at a = 0, hence the 
absence of singlet-SC order. However, the triplet states 
have nontrivial minima. It has been argued[12] that, in 
the presence of disconnected Fermi surfaces on the trian- 
gle or honeycomb lattices, triplet pairing may be favored 
over singlet pairing. 

The nodal structure of each of the four solutions (7- 
10) is illustrated on Fig. 3 in plots of the order parameter 
{cA,ii,TCA,-k,i), taken from corresponding mean-field so- 
lutions. Solution (7), which is a pure Bi representation, 
has three nodal lines around each Dirac point and may 
be thus qualified as f-wave. The other three all have a Ei 
component and display a single nodal line: they are p- 
wave. The last solution (10) is complex and only its real 
part is represented; its imaginary part has a similar struc- 
ture with conjugated nodes and its modulus is constant 
around each Dirac point; it is a p-t-ip solution. It can 
be shown that if only NN pairing is present in the mean- 
field solution of type (7) , then the dispersion relation (3) 
amounts to a simple rescaling of the noninteracting case 
(a renormalization of the Fermi velocity) and supercon- 
ductivity is a hidden order. This is due to the pairing 
function ?7k being proportional to the hopping function 
7k in that case and is an artefact of the restriction to NN 
pairing. In reality, pairing extends to further neighbors 
(even if NN pairing only is used as a Weiss field in VGA) . 
For this reason, we added a third-neighbor pairing term 
to the mean-field Hamiltonian of solution (7) in order to 
produce the top-left plot on Fig. 3. 

For each value of the chemical potential /it, a normal- 
state solution can also be found, by suppressing the Weiss 
fields. The values of the Potthoff functional for these 
different solutions gives an estimate of the energy density 
E — + fj,n, as a function of electron density n = Tr G. 
The top panels of Fig. 4 show the condensation energy 
(i.e. the energy of the normal solution, minus that of the 
ordered solution) for each of the four SC solutions as a 
function of doping S. We see that the p-wave solutions 
(8), (9) and (10) are close in energy to each other and 
seem to be favored at larger coupling (except close to 
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FIG. 4: (Color online) Top panels : Condensation energy (in 
units of t) for the various triplet superconducting solutions as 
a function of hole doping {5 = 1 — n), obtained with 6-site 
(left) and 10-site (right) clusters, at U — 6. Bottom panels: 
Root-mean square order parameter for the same sublattice 
(left) and different sublattices (AB), as a function of doping. 



half-filling). As an additional measure of the strength 
of the SC solutions, we computed the root-mean square 
(RMS) value A(:'j„s of the momentum-dependent order 
parameter: 
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where fi and v are sublattice indices {A or B). This 
is illustrated on the bottom panels of Fig. 4 for the 6- 
site cluster. Note that the order parameter falls towards 
half-filling, and that the diagonal piece falls to zero. 
The diagonal and off-diagonal parts are similar, except 
for the f-wave solution (7). A condensation energy of 
O.OU corresponds, in this case, to an energy scale roughly 
equal to room temperature, but infinite-range order in 
a two-dimensional system would not exist at a nonzero 
temperature. 

The possibility of ferromagnetic Neel order, arising 
from the bipartite character of the lattice, was also inves- 
tigated (this order is defined at zero wavevector (Q — 0), 
hence the term 'ferromagnetic'). A Weiss field M, mul- 
tiplying the staggered magnetization operator 
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was also treated as a variational parameter. No such or- 
der was found away from half-filling. At half-filling, a 
Neel order was found for U > 3, and that solution is en- 
ergetically favored over the SC solutions for U > 6.5. As 



said above, the relevance of the (non-extended) Hubbard 
model at half-filling in this sytem is unclear. However, 
away from half-filling, the exchange of ferromagnetic fluc- 
tuations is a possible mechanism of triplet pairing. 

Thus, our second conclusion is that triplet supercon- 
ductivity exists in the doped system, and that a p-wave 
solution seems favored. It is impossible to reliably state 
which of the three solutions (8,9,10) is preferred in the 
thermodynamic limit. In the case where solution (10) 
is preferred, time-reversal invariance would be sponta- 
neously broken. 

Remains the issue of the physical realization of 
graphene sheets with sufficient doping for this prediction 
to be tested. At present, only very small doping can be 
achieved by applying electric fields {S ^ lO^''). DFT cal- 
culations show that coating graphene with a metal can 
induce a shift in chemical potential of '--^ 0.5eV [13], which 
translates, in the non-interacting case, into S ~ 0.03, but 
which may translate into an even smaller value if corre- 
lations are taken into account. Chemical doping of some 
kind may be the only possible way to reach values of 6 
relevant to this work. 
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